Choosing the correct frequencies

import numpy as np
import corner
import pandas as pd
import matplotlib.pyplot as plt
import exoplanet as xo
import pymc3 as pm

%config IPython.matplotlib.backend = "retina"
from matplotlib import rcParams
rcParams["figure.dpi"] = 100
rcParams["savefig.dpi"] = 150

The default behaviour in Maelstrom is to estimate a set of frequencies for the model if none are supplied. While in most cases, these frequencies are fine, there are cases where special care must be taken. Let’s first generate a synthetic light curve with TESS sampling. We’ll consider a very rapidly oscillating star (frequency of 300 cpd), in a binary system with period 5 d. This is a rather unusual system, but we’ll use it to prove a point.

from maelstrom.synthetic import SyntheticBinary

period = 5.
eccen = 0.5
asini = 100
varpi = 1.
tref = 0.

n_freqs = 1
freqs = np.array([300.])
amps = np.array([1])

time = np.arange(0, 27, 1.0 / (24 * 30))

lc = SyntheticBinary(time, freqs,amps,
                     period, eccen, asini, varpi, tref)

lc.add_noise(snr=500)

Now, let’s ask Maelstrom to have a quick look at the star for any time delay signal.

from maelstrom import Maelstrom
ms = Maelstrom(lc.time, lc.flux / 1e3,
              max_peaks=1
              )
ms.first_look(segment_size=0.1);

What’s going on here? You’ll notice that the time delay signal seems to be steadily decreasing over time. We can see the binary motion, but if we try and fit it we’ll get some nonsense back out. Maelstrom wasn’t designed for this!

In fact, this motion is caused by Maelstrom incorrectly picking out the correct frequency. Frequency modulation causes a series of sidelobes to manifest in the amplitude spectrum. In some cases, such as stars with high frequency and low orbital period, these sidelobes can have higher amplitudes than the true frequency. So, let’s take a closer look at the synthetic data by calculating the amplitude spectrum.

fmin, fmax = freqs[0]-3*(1/period), freqs[0]+3*(1/period)
ax = ms.plot_periodogram()
plt.xlim(fmin, fmax)

ax.axvline(ms.freq, c='r', linestyle='dashed');

You’ll notice that the true frequency at 300 cpd is split on both sides by a series of peaks. These peaks separated by the orbital frequency of the binary system (1/Porb). When Maelstrom looks for the frequency, it grabs the highest peak and optimizes that. Indeed, if we grab the sidelobe immediately to the right of the central peak, we’ll notice that the time delay signal is instead increasing over time:

from maelstrom import Maelstrom
ms = Maelstrom(lc.time, lc.flux,
              freq=[300.2]
              )
ms.first_look(segment_size=0.05);

Hopefully this illustrates how important it is to choose your oscillation modes carefully. Finally, let’s make sure we can recover the original time delay signal by picking the correct frequency:

from maelstrom import Maelstrom
ms = Maelstrom(lc.time, lc.flux,
               freq=freqs,
              )
ms.first_look(segment_size=0.05);